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1-C Abstract 

Q Large scale structure surveys are likely the next leading probe of cosniological information. It is there- 

-M fore crucial to reliably predict their observables. The Effective Field Theory of Large Scale Structures 

j^ (EFTofLSS) provides a manifestly convergent perturbation theory for the weakly non-linear regime, 

where dark matter correlation functions are computed in an expansion of the wavenumber k over the 

wavenumber associated to the non-linear scale /cnl- To push the predictions to higher wavenumbers, it 

is necessary to compute the 2-loop matter power spectrum. For equal-time correlators, exactly as with 



> 






■^ standard perturturbation theory, there are IR divergences present in each diagram that cancel completely 



in the final result. We develop a method by which all 2-loop diagrams are computed as one integral, 
with an integrand that is manifestly free of any IR divergences. This allows us to compute the 2-loop 



(^ power spectra in a reliable way that is much less numerically challenging than standard techniques. We 

^^ apply our method to scaling universes where the linear power spectrum is a single power law of k, and 

• • where IR divergences can particularly easily interfere with accurate evaluation of loop corrections if not 

. !-H handled carefully. We show that our results are independent of IR cutoff and, after renormalization, of 

^\ the UV cutoff, and comment how the method presented here naturally generalizes to higher loops. 

S 



1 Introduction 

Large scale structures have the potential of becoming the next leading cosmological probe for the 
initial conditions of the universe. Our knowledge of the initial conditions scales as the cube of the 
maximum wavenumber in observations that we can theoretically predict. This is a tremendous 
amount of potential data. 

On short scales density perturbations have become non-linear and have undergone collapse. 
A description of this regime likely necessitates the use of A^-body numerical simulations, with 
analytical techniques providing guidance. The situation is very different on large scales, where 
non-linearities are weak and an analytical treatment should be possible. The recently formulated 
Effective Field Theory of Large Scale Structures (EFTofLSS) [1, 2] is the theory that allows 
us to consistently make predictions for correlation functions at a certain wavenumber k in an 
expansion in 6{k) = 6p{k)/p or equivalently in powers of k/k^i^, with /cnl being the wavenumber 
corresponding to the non-linear scale. Since we are interested in computing correlation functions 
for 6{k) <C 1 or equivalently k ^ A;nl, this perturbation theory is manifestly convergent. The 
effective theory differs from standard, normally used, perturbation theories [3] by additional terms 
in the fluid equations of motion, such as the speed of sound, viscosity, and stochastic pressure. In 
addition, loops are performed formally with a cutoff and these additional terms of the EFT have 
the role of canceling the cutoff dependence for physical observables and to add a finite contribution. 
After this step, the theory is said to be renormalized, and it is only at this point that the expansion 
in fc/^NL is manifest [2]. The EFTofLSS was developed and applied to simulation data in [2], where 
it was explicitly shown how the additional terms of the EFT remove the cutoff dependence of the 
loops and how the theory is renormalized. Since the additional terms such as the speed of sound 
are incalculable within the effective theory, they need to be either measured from observed data, 
or extracted from A^-body simulations. Both approaches were developed in [2] , where it was shown 
that the 1-loop prediction of the EFTofLSS is in agreement with the dark matter power spectrum 
to within a percent up to the relatively high scale of A; ~ 0.24 h Mpc~^ . The remarkable agreement 
with observations, the improvement with respect to other currently available techniques (see for 
example [4]), and the self-consistency of the different methods of extracting the parameters of 
the EFT, give very strong evidence that the EFT is the correct language to make theoretical 
predictions for LSS. This becomes even more evident when one tries to make predictions for other 
toy, 'scaling', universes where the initial power spectrum follows a simple power law. In this 
set up, all techniques other than the EFTofLSS are unable to make predictions, while the EFT 
can [5, 6]. By using the renormalization techniques developed in [2], Ref. [6] has explicitly verified 
that indeed the EFTofLSS is able to make predictions for the scaling universes. 

From the perspective of the EFTofLSS, the former techniques are missing important terms 
in the equations of motion. Once these terms are included, perturbation theory may converge 
for much larger k than previously thought. Given the importance of increasing the window of 
modes over which we can reliably predict observables, and given the encouraging results of [2], 
it is tempting to perform calculations beyond one-loop. While conceptually straightforward, the 
next-to-leading order contribution to the power spectrum, from 2-loop diagrams, happens to be 



technically subtle. Apart for the presence of counterterms that need to be evaluated at lower loop 
order, the calculation of the 2-loop contribution is actually exactly the same as in Standard Per- 
turbation Theory (SPT) [3]. The calculation is numerically challenging. There are four diagrams, 
denoted as P51, P42, P33 and P33 , each one involving a 5-dimensional integral. One of the main 
difficulties is that each diagram has infrared (IR) divergences as some of internal momenta go to 
zero. These must cancel between diagrams given infinite precision, but finite resources makes it 
easy to generate spurious numerical errors. It is a peculiar fact associated to the nice IR behav- 
ior of the standard real-universe linear-power spectra that such spurious issues have largely gone 
unnoticed until now - although as we will discuss they are very much a real issue with noticeable 
consequences. These issues become more evident in the case of scaling universes that we treat 
here. 

As proven in [7, 8], in equal-time matter correlation functions, all IR divergences must cancel 
when we sum together all contributing diagrams, provided that the linear matter power spec- 
trum Pii grows in the IR more slowly than fc" with n > —3 (see [9] for a recent discussion). 
The physical reason behind this is the following. If we consider two short modes, they will see 
long modes as a background, and they will be carried over by the velocity of the long modes, the 
so-called bulk fiows. By the equivalence principle, the dynamics of the short modes is sensitive to 
long modes only through second derivatives of the metric, the so-called tidal forces, which go as 
k^Pii{k). If 72 < —3, these tidal forces grow larger for longer wavelengths, leading to IR divergent 
correlation functions. For equal-time correlation functions this is all that there is to it, and so 
IR divergences are guaranteed to be absent for n > —3 by this argument. For non-equal times 
correlation functions, there is an additional effect through which IR modes appear. This is in 
inducing a relative displacement between the fluctuation evaluated at the earlier time and the one 
at the later time. Since the power spectrum of velocities goes as k^Pii{k)/k'^, these IR divergences 
at non-equal times will cancel only for n > —1. 

In this paper we will focus on equal time correlations with n > —3, where IR divergences are 
guaranteed to cancel. However, the fact that this cancellation happens only when summing over all 
the diagrams introduces difficulties at the numerical level. In fact, in the currently available codes 
that perform the two-loop calculation for SPT and that we have tested [4, 10], each diagram is 
evaluated separately. This means that each diagram is IR divergent and so needs to be regularized 
in the IR. Ideally, if the integrals are computed exactly, the large IR dependent part will cancel and 
we will be left with a smaller, IR- in dependent part. However, this means that the integrals need 
to be evaluated to a much greater precision than what is really necessary, because the leading part 
will actually cancel. To be explicit, if we wish to compute the two-loop power spectrum with 10~^ 
relative precision, and the IR dependent part of each diagram is 10^ larger than the finite part, 
we need to evaluate each diagram to 10"'' precision. This is clearly a numerical challenge and a 
waste of resources. Even worse, as we will explain, in an expansion of the IR-independent part 
in powers of /c/Zcnl, the leading IR-independent part is actually degenerate with the effect of a 
counterterm from the EFTofLSS, the speed of sound term. This means that we are not interested 
in computing the leading IR-independent part, but actually the next subleading part. This raises 
the need for higher, numerical precision, making the challenge even harder. 



We find that the currently available 2-loop codes that we have tested [4, 10] are not precise 
enough for scaling universes, probably for the aforementioned reasons ^. In contrast we present 
and implement what we call the IR-safe global integrand. The main idea will be the following: 
since IR divergences cancel between diagrams, we will bring the integrands of all the four diagrams 
under one integral. Additionally, in order to cancel IR divergences that happen in different regions 
of the phase space (for example as the internal momentum goes to zero or to minus the external 
momentum), we will perform some changes of variables in the integration coordinates so that all 
the IR divergences happen only as the internal momenta go to zero. In this way, the resulting 
integrand will be IR convergent, the cancellation of IR-divergences will not rely on numerical 
precision, but it will happen directly at the level of the integrand and therefore will be automatic. 
Consequently, the integration computes the IR-finite part directly, so that the precision with which 
the integrals are evaluated is the precision of the IR-independent part. 

Because of the numerical challenge, we will apply and test our approach on a scaling universe. 
This is interesting and useful also because simple dimensional analysis allows us to estimate the 
parametric result of the integral, allowing a direct verification of success. Furthermore, scaling 
universes do not have a property of the current universe that gives the illusion that a code, when 
applied on the true universe, is IR-safe. In fact, in the current universe, for scales longer than the 
equality scale k^q, the power spectrum scales as k^, which makes each loop diagrams much less 
IR divergent. Because of this, it is possible that some numerical approaches only approximately 
cancel the IR divergences, but do so well enough for scales longer than the equality scale, so that 
most of the contribution for a given high-fc mode comes from keq- This result would be naively 
IR-independent, but it is still not the correct answer, as, for correlations evaluated at equal times, 
every k modes receives contribution only from its neighboring or higher-/c modes. In a separate 
paper [11], we will apply the IR- safe global integrand to give the 2.5-loop prediction of the power 
spectrum from the EFTofLSS, where we will see that the EFTofLSS predicts the power spectrum 
with percent precision up to roughly k ~ 0.5 hMpc~^. 

2 The IR-safe global integrand 

In this section we will be interested only in the loop terms, so we will therefore momentarily 
neglect the EFT counterterms. This means that the expressions we write in this section are the 
same as in SPT. Furthermore, when necessary to use explicit linear power spectra, we will consider 
scaling universes where the linear power spectrum, defined as P = Pn, is taken to be 



^«^i(£i • « 



with the tilt, n, being some number. 



"'^We stress the codes we will use here for comparison are not designed for scaling universes. It might well be 
that some of these IR-issues are numerically irrelevant for these codes in the case of the true universe. 



2.1 1-loop Integrand 

The 1-loop matter power spectrum is often written as 
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-Pl3 + P22 
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with the kernels F2 {q,k — q) and F3 (k, q, —q) given by 
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As g — )■ 0, the integrands have the following asymptotic behavior (ignoring a common factor of 
(27r)"^ and including a g^ factor from the integration in spherical coordinates): 



q'pi3ik,q) ~ -P+VV + C?(g 



n+2\ 



q'^P22{k,q) 



I 

g-^O 2 



k^+^H^q^ + C(g 



,n+l^ 



(5) 



where n = {k ■ q)/{kq). We will discuss the 0{q"'~^^) term momentarily. One sees that the leading 
IR divergence in P22 does not cancel the one in P13. In fact, there is another IR divergence in P22 
as g*— )■ A;, associated with sending the internal momentum q — k to zero: 



g'P22{k,q-) ~4^"^^rfm^i^-^T + c?(i^-r^^: 

q~^k ^ \k — q\ 



(6) 
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Remarkably, if we sum the leading IR-divergent contributions, they cancel and we are left with an 
IR independent result. This is guaranteed to happen from the Galilean invariance of the equations 
of motion [7, 8], and is also guaranteed for all sub leading IR divergences. However, the presence 
of the IR divergences in each of the diagrams makes the calculation numerically much more 
challenging. If one is to evaluate P13 and P22 separately and then sum them up at the end, one 
needs to put an IR cutoff, compute two very large numbers, and then add the partial results only 
to obtain a subleading piece. This affects the precision of the result: a certain relative precision 
in Pi3 or in P22 translates into a much worse relative precision in the final result for Pi_ioop- 

One might think that this problem appears only in the toy scaling universe that we are consid- 
ering. In the true universe, Pn ~ k for k < keq, making the IR divergence of P13 and P22 vanish. 
Still, this means that the result of P13 and P22 is dominated by the contribution of modes of order 
keq, whose contribution will cancel for modes k ^ keq- While at 1-loop the numerical challenge is 
relative, and one can afford such an inefficiency, the challenge becomes much harder for the 2-loop 
case that we consider next. 

Since IR divergences are guaranteed to cancel, it would be ideal to compute the 1-loop (or 
any n-loops) power spectrum without IR divergences ever appearing. Ideally, one would like to 
have an integrand that does not have any IR divergence. Clearly, given that the IR cancellation 
happens between P22 and P13, the only possibility is to write the two as an integration of a single 
integrand: 

Pl-loop IR-safe guess = / 7^"^ \Pl^{k,q) +P22{k,q)\^ ■ (7) 

This is not enough though, as P22 diverges in a symmetric way, both as g — ?■ and as g — )■ A;. 
However, it is possible to perform the following change of variables on P22 in order for it not to 
have any divergence aX q ^ k. The idea is split the region of integration in order to isolate each 



singular term, and to map the singularity at g—)- A; to one at g — )■ 0. We write 



P, 



22 



d?q 



P22{k,q) 



d?q 

\q\>\k-q\ (271" 

d?q 



3 P22{k,q) 
P22{k,k-q) 



(9) 



.^ .^P22{k,q)+ , 

\<^<\k-q\ {^T^) J\q\<\k-i\ (2vr)^ 

= 2/ . ^^P22{k,q)=2 [^^p,,{k,q)Q{\k-q]-q) 

J\q\<\k-g\ {^V J {^^> 

where in the first passage we have split the region of integration, in the second we have changed 
integration coordinates from q to q = k — q, and in the third we have relabelled q as q and used 
the property P22{k, q) = P22{k, k — q), obtaining a factor of 2. 

In this form, P22 has an IR divergence only as g* — ;■ 0, and we can see from Eq. (5) that it 
will cancel with the corresponding term in pi^. The final important point is that P22 also has a 
subleading divergence (the 0{q"'^^) piece in Eq. (5)), with the following behavior: 



q^P22{k,q) ~ ^A;"+VV + A^"^V(6 + 8/i2-7V)g"+i + 0(g"+2) 
<?-s>o 2 14 



(10) 



This term, being odd in /i, disappears once the angular integral is performed, and can also be 
eliminated in the integrand by symmetrizing in q and —q. Thus, we can write the final IR-safe 
-^ i-ioop as 



Pi- 



loop IR-safc 



d^q 
(2^ 

d^q 
(27r)3 



Pi-i{k, q) + P22{k, q)Q{\k - q] - q) + P22{k, ~q)Q{\k + q\ - q) 

6P{k)P{q)F^'\k,q,-q) 

+2P{q)P{\k - q\) [Ff (g, k - q}]' e{\k - q\ - q) 

+2P(g)P(|fc + g-j) k^^)(-g-', fc + g-)!' e(|fc + gl - g) 



^One can derive this more simply by un-doing the Dirac (5- function of momentum conservation: 
P22 = 



—^ P22ik,q) = / 7^ / tSt Sn{q + P + k) 2P{q)P{p) \F^'\q,p)Y 



(27r)3 ^^^^■■^'■^ J {2tt)3 J (27r)3 
d^q f (Pp 



(11) 



(8) 



{2ttY J (27r)3 
cPq 



5T,{q + p + k) Q{p~q)2P{q)P{p) F^'> {q,p) 



?i^)i 



(27r)3 



e(|fc - gl - q)2P{q)P{\k - q\) F^ {q, k - q) 



,i-)l 



2 / ^P22{k,q)e{\k~q\~q)., 



where in the first passage of the second line we have used the symmetry of the integrand under exchange of q and p. 



There are two advantages to writing the 1-loop integral in this way. A trivial one is that we 
have to do one integration instead of two, which is less time-consuming. More importantly, a 
second one is that the integrand has no divergences as g — ?■ for any n > —3. This is guaranteed 
to be so: the sum of all IR divergences must cancel in the integral, and since we have only one 
integrand and the internal momenta go to zero only for g*— )■ 0, the integrand itself must be finite, 
or at least integrable, for n > —3. As we have just explained, this makes an accurate numerical 
integration much easier to perform. These manipulations of the integrand are not so important for 
the (relatively straightforward) 1-loop integrals, but become essential when one tries to evaluate 
the much more challenging 2-loop diagrams. 



2.2 2-loop Integrand 

Now, let us examine the 2-loop integrals. They are usually written as 
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where the kernels -^345 can be calculated from the recurrence relations found, e.g., in [3], and 
symmetrizing over all arguments. We repeat these relations here for convenience: 
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where fci = gl + ■ ■ ■ + g^, ^2 = Qm+i + ■ ■ ■ + Qn, k = ki + k2, and F„ = (?„ = 1. From now on, we 
define the integrands psi, P42, P33 by 

^''^^^ ^ J (2^ J ~(2^ ^''^^' ^"'^ ' ^^^^ 

and we also use the notation jj, = [k ■ q)/{kq) and u = {k ■ 'p)/{kp). 

We want to write all the 2-loop integrals as only one integral where the integrand is a sum 
of terms, each one with its leading divergence located only atg*— J-O&p— )-0. In addition, 
each diagram will have subleading divergences occurring when only one combination of internal 
momenta tends to zero (for example, g — )■ & p fixed, or p — )■ & g fixed), and we would like 
to map all such divergences to g — ?■ & p fixed. Here, we examine each of the four integrands 
separately: 

1. P51 only has a single leading IR divergence, at g — )■ & p — )■ 0, so it already satisfies our 
first condition. It also has subleading divergences when either q or p approach zero while 
the other one remains fixed. Since P5i{k,q,p) is symmetric in g'and p, we can perform the 
following manipulations to move the divergence at p — ;■ to g — )■ 0: 

It is also possible to find the explicit forms of the various divergences. The leading divergence 
takes the following form (omitting a factor of (27r)~^ and including the factor g^p^ from the 
integration measure, since it affects the degree of divergence of the integral): 

gV2p5i(fc,g-',p1. ~ lk'+y'uYp\ (17) 

while both the leading and subleading behavior is captured^ in the following expression: 

qp2p,,{k,q,p)^^^ 42([p+p2p_4^2p2^2) ^P ' 

(18) 
with further subleading divergences discussed later. We will find that both the leading and 
subleading divergences will cancel when summed with the other diagrams. 



•^Deriving this expression analytically is not straightforward given the large number of terms in psi. One can 
instead check that they agree numerically. 



2. The leading divergences of P42 are atp— J-O&g— J-O, and p— i-O&g— i-A;, sowe must 
manipulate the integrand to re-map the latter divergence top— J-O&g— i-O. The algebra 
works the same as for P22: split the q integral into q < \k — q\ and q > \k — q\ pieces, and 
use the substitution q = k — q on the second piece. We end up with 

^^'^^^ = J^sJ^3^P^^(^^'l^P)^(\^-'^-l)^ (19) 

which, however, has a subleading divergence at p — )■ in addition to the one at g* — )■ 0. 
To remedy this, we first symmetrize the integrand in q and p (recall from Eq. (13) that 
PA2{k,q,p) is not symmetric in q and p), then perform the same steps as in Eq. (23) to 
restrict the domain of integration to p > g. The result is 

P42{k) = / (03 / (03 2 [P^2{k, q,p) Q{\k -q]-q) + {q ^ p)] Q{p - q) . (20) 

The leading divergence is 

q^pHp42{k,q,p) ~ -2A;^+>VVp" • (21) 

The subleading divergences for p42{k, q, p) and p^iik, p, q) must be treated separately because 
they contribute over different domains in p. One finds 

gV2p42(fc,g,p) ~ 

A;^+V(21A;V2 +p4[_io + 59^/2 _ 2%v^] - 2A;V[5 - 22//^ + 38z/^]) „ „ 

21([A;2+p2]2_4Pp2^2) ^ P 

22 /T ^ ^ A;6+V(7A;i^+p[3-10i/2])2 



>(I) 



3. P33 has its leading divergence at p — )■ & g — )■ 0, and subleading divergences at g — )• and 
p — >■ 0. The can therefore be rewritten in the same way as P51: 

Indeed, it diverges in the same way as P51, with leading piece 

gV2pg(^,g-',p)^ ~ V>VgV. (24) 

ij-!>0,p-s>0 Z 

and subleading expression 

gV 2py(^,g,p) r.^ 42([P+p2]2_4,U2) ^^> 

(25) 
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4. P33 ^^^ three leading IR divergences, atg— T-O&p— 7>0, g— 7>A;&p— t-O and g — )■ & 
p ^ k. There are also three subleading divergences: g*— )■ & p fixed, p — > & g* fixed, and 
•p -^ k — q. Thankfully, these can all be handled in a systematic way. First, we rewrite the 
integral to be symmetric in g*, p, and a third internal momentum t. 

xPuiq)Piiip)Piii\k-q-p\) 
<Pp f d^q f d^£ 



j-^M^\imF^\-^.-p. 



(27r)3 J (27r)3 J (27r)3 

xPn(g)Pii(p)PiiW5D(^-g"'-p-/) . (26) 

Written this way, the leading divergences occur when two of g*, p, and d, go to zero, and the 
subleading divergences occur when only one momentum does. We can split the region of 
integration using the following sum of step functions: 

0(p - g) [0(£ -p) + e{p - £)e{q -£) + e{p - £)&{£ -q)] + {q^p) . (27) 

Since the integrand is symmetric in all three momenta, we are free to relabel the momenta 
in each term; after doing so, we find that each one is equivalent, so that we can write the 
sum as 60 (p — g)0(^ — p), and the integral as 

xP^,{q)Pu{p)Pnii)SB{k - q- p- £) x6Q{p- q)Q{£ - p) . (28) 

This ensures that the leading and subleading divergences are all mapped to the same loca- 
tions: (g— )-0&p— )-Oorg— )-0& p fixed), as desired. Upon evaluating the delta function, 
we are left with the original integrand P33 times a pair of step functions: 

Written like this, the integrand's leading IR behavior is 

gV6p£)(^,g-',p-)^ ~ A;^+VA>", (30) 

while the subleading behavior is 

6g P Pss {Kq,pl -^ 49(P+p^-2A;H^ ^ ^ " ^''^ 
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In summary, we have rewritten the 2-loop integrand in the following way: 

P2-ioop{k, g, p) = 1 2p5i (fc, q,p) + 2p^5 (^> ^, P) 

+2 pi2ik,q,p)Q{\k-q\ - q) + P42ik,p,q)Q{\k - p\ - p) 
+6pg\k,q,p)Q{\k - q- p\ - p)^ Q{p - q) , (32) 

or, rewritten in terms of SPT kernels and linear power spectra, 

P2-iooAk,q,P} = {60F^'\k,q,-q,p,-p}Pnik)Pniq)Pnip) 

+18F^''\k,q,-q)F^'\-k,p,-p)Pu{k)Pu{q)Pu{p) 
+48Fi'\q,k-q)Ft\-q,q- k,p, -p)P,,{q)P,,{p)P,,{\k - q]) Q{\k - q] - q) 

+A8Fi'\p,k - p)Fi'\-p,p- k,q,-q}Pniq)Puip)Pii{\k - p\)ei\k - p\ - p) 

+3QF^''\q,p, k-q- p)F^^\-q, -p, -k + q + p) 

^P^MPii{p)Piii\k-q-p\)Q{\k-q-p\-p)]Q{p-q) . (33) 

By summing Eqs. (17), (21), (24), and (30), we see that the leading IR divergences have cancelled. 
In addition, by examining Eqs. (18), (22), (25), and (31), we see that the subleading divergences 
with asymptotic behavior g" as g* — )■ & p fixed have also cancelled. In particular, we can notice 
that in the limit oi q <t^ k, p, the structure of the 9-functions forces two subsets of the contributions 
to cancel independently among themselves. In practice, we have 

P2-ioop{k,q,p) ~ <S>{p-q) X (34) 

{P^^{k)P^r{q)P^M X 

(eOFf (fc, g, -g-p, -p) + 18Ff (fc, q, -q)F^'\-kJ, -p) + 48Fi^)(g-', k - q)Ft\-q, q- k,p, -p)) 

+Q{\k-p\-p) Pu{q)Puip)Piii\k-p\) X 

(48Fi^^ (p, k - p)Ff\-p, p-k, q, -q) + 36^^^^ (g"*, p, ^ - g"*- p)F^\-q, -p,-k + q + p))] . 

Since for |fc — p| — p < the terms in the fourth and fifth line vanish, this means that the IR 
divergences from the terms in the third lines cancel independently among themselves. This in 
turn implies that the terms in the fifth line cancel independently among themselves. 

For n > —2, there are no remaining IR divergences, since terms scaling like g""*"^ or higher 
will converge when integrated. For n < —2, terms with ~g"+^ asymptotic behavior will result in 
divergences. However, these terms all involve odd powers of /i or i/ (recall that these represent 
k ■ q and k -p respectively), which vanish when the angular integrals are performed. Alternatively, 
these terms vanish when the integrand is symmetrized in g'& — ifand p & —q, so this is the final 
step we must take to ensure cancellation of divergences in the integrand itself for n > —3. 
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Therefore, the following form of P2-I00P is free of IR divergences before integration: 

/d v f do 1 r -> -> 

J2tt)^ / f27r)3 4 r2-ioop(/^, q,P)+ p2-\oop{k, -g, p) 

+P2-ioop{k, q, -p) + p2-ioop{k, -q, -p) , (35) 

where P2-\oop{k,q,p) is given in Eqs. (32) and (33). 

As g* — )■ & p — )■ 0, the remaining integrand scales with momenta like g"-+2p"'+2^ while for 
g — )■ with p fixed, it scales like g""*"^. As previously discussed, this makes an accurate numerical 
calculation of P2-I00P much easier, as the delicate cancellation of large values of the integrand at 
different locations in momentum space is no longer necessary. 

2.3 Higher Loops 

IR-divergencies will cancel at any loop order for equal time correlators [8]. It is to be expected 
that the procedure we have outlined at 1-loop and 2-loops to define a global IR-safe integrand 
should work relatively unchanged at higher loops. The main steps of the procedure should be the 
following. As naively written, inside L-loop integrands the linear power-spectra will be functions 
of linear combinations of L-loop momenta and external momenta. The first step is to explicitly 
introduce one momenta-conserving Dirac ^^-function so that all linear power-spectra that appear 
in the integrand are functions only of individually labeled momenta {g^j's, and not of combi- 
nations of (ji among themselves or with k *. Then, one symmetrizes the integrand with respect 
to permutations of the subset {g!,} of the {gi} that are involved in the Dirac (5D-function, and 
also with respect to all parity transformations acting on the various {g^}'s. By inserting proper 
G-functions and symmetry factors, one g-orders the various integrations over the moduli of the 
{g^}'s (similarly to the time-ordering procedure familiar in Quantum Field Theory). In this way 
one ensures that only one of the {g^} momenta is allowed to go to zero. At this point, one sym- 
metrizes with respect to all permutations including the remaining {gj}'s and all parities acting on 
those, and orders the resulting symmetric expression with B-functions and symmetry factors so 
that only one momenta out of the full set of the {gjj's is allowed to go to zero. As the last step, 
one evaluates the introduced Dirac ^^^-function, so that the dimensionality of the final integral is 
no worse than when one started. 

''why only one additional Dirac (J^j-function necessary no-matter how high you go in loop order? Because each 
diagram is an L + 1 primordial-fluctuations sewing of two tree-diagrams, so only one "cut" momenta ever need be 
written as a linear combination of the others, and it is that one that we resolve by introducing a Dirac (Ju-function. 
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3 Application to Scaling Universes 

3.1 Power Counting and Perturbative Expansion 

Let us apply our IR-safe integrand to the computation of the matter power spectrum at 2-loops 
in the case of a scahng universe where the tree level power spectrum 

We choose the scaling universes because simple dimensional analysis allows us to infer the k- 
dependence of the answer. We will specialise to the cases n = —3/2 and n = —5/2 when doing 
the numerical calculations. 

The first point to notice is that loop integrals will diverge in the ultraviolet (UV), the details 
depending on the tilt n. This simply means that the divergent contribution needs to be cancelled 
by counterterms that make the answer finite and cutoff-independent. The EFTofLSS provides the 
counteterms that cancel all physically possible UV divergences. This was explained and applied 
to the current ACDM universe at 1-loop in [2] , and later at 1-loop for the scaling universes in [6] . 

Following [12], if L is the number of loops, I the number of internal lines, V the number of 
vertexes touching loops, the superficial degree of divergence D of a diagram is given, in terms of 
the tilt n of the linear power spectrum, by 

D<3L + nI-V . (37) 

The factor of 3L comes from the integration measure, the factor of nl comes from the momentum- 
scaling of the internal lines, while the factor of — l^ comes from the UV-properties of the integrand 
F^^^ that vanishes at least as 1/g, with q being one of the high-scale internal momenta. This 
decoupling property is stronger when only one of the internal momenta goes to infinity. In that 
case the contribution from a vertex goes as 1/g^, but this property does not hold generically when 
more than one momentum goes to infinity. 

Let us call Dg the superficial degree of divergence of the loop diagram Pg, where the label g 
is given by a couple of numbers, each one representing the perturbative order of two contributing 
(5's in the diagram. In order to cancel the IR divergences, we need to consider all diagrams at 
the same time, but the diagrammatic splitting is useful for studying the UV part, as this dictates 
which kind of counterterm of the EFTofLSS will remove the divergence. At 1-loop we have 

Di3 = n + 1 , D22 = 2n + 1 (38) 

We see that P13 can diverge forn > —1, while P22 does so for steeper power spectra with n > —1/2. 
At two loops, we have 

DlS = 272 + 5 , ^24,33 = 372 + 4 . (39) 

We see that P15 can diverge for n > —5/2, while P24 and P33 will generically diverge for steeper 
power spectra with n > —4/3. 
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The superficial degree of divergence represents an upper bound to tlie degree of divergence of 
a diagram. For example, it does not even take into account constraints from rotational invariance. 
Even more powerfully, one can check if a potential divergence can or cannot be reabsorbed by one 
of the counterterms of the EFTofLSS to further constrain the structure of the divergent terms. 
For example, for n > —5/2 and for P51, one can realize that the dependence on /cnl is fixed to be 
^NL " ; then use the superficial degree of divergence and finally the fact that P51 has one Pn out 
of the loop-integral, to conclude that P51, evaluated with a UV cutoff A, takes the following form 

P51 = (40) 

where q are order one numbers. If this answer were correct, there would need to be a local coun- 
terterm in the EFTofLSS that goes as /c x Pn ~ k{6l). Since the power spectrum is proportional 
to Pii(/c), the counterterm must be proportional to 6. Such a counterterm should therefore have 
the form of k 6, but this breaks rotational invariance or locality (which requires analyticity in k) 
and so is absent from the EFTofLSS. The lowest derivative available counterterm is the so-called 
speed of sound term of the form k'^6, which indeed can reabsorb the subleading divergence ^. We 
can therefore drop the most divergent term, and conclude that 



Con+A I 1 — I ( 1 — I Pn + subleading divergent terms . (42) 



For n = —3/2, the 1-loop terms are finite, but the 2-loop terms can (and actually will) diverge 
as A^. In addition, there could be in principle a subleading logarithmic divergence of the form 



A\ / k ^^ 



n.„=-3/.3cf„,log^^j ^-1 Pu, (43) 

for which there is no available counterterm after using the rotational invariance of the counterterm. 
We therefore conclude that the logarithmic divergence is absent: c^ = 0. One can arrive at the 
same conclusion starting directly from the loop integral by expanding in fc/gioop, where gioop — ?■ A 
is any of the diverging loop momenta. The first sub-leading term (after the leading divergence) 
must be of the form /(gioop)^ioop ■ ^/o'loop which vanishes after performing the angular intergrals. 
We see that the EFTofLSS allows us to infer information about the outcome of the integral even 
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The same identical reasoning would have allowed us to remove a putative divergent term of the form 

A \ 2n+6 
A / A 



without having to deal with the UV properties of the vertex functions. Indeed, this term can be trivially excluded 
due to the absence of a proper counterterm in the EFTofLSS. 
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before actually doing the calculation and beyond naive dimensional analysis. It is interesting to 
verify these claims in a numerical calculation as we do next. 
Although the logarithmic term is forbidden, a finite term 



k ^' 



P5I, n=-3/2 D Cfinite 1 ( ^ ) ^H (44) 

is allowed with Cfinite i 7^ 0, since it cannot forbidden on the grounds of available counter-terms. 
Because the divergences will all be absorbed by the counterterms, the finite terms are the only 
physically meaningful contribution to the power spectrum. In this case there is no scale in the 
computation except for the external wavenumber k (the overall dependence on fc^L is irrelevant as 
it simply enter linearly in the answer), for the loop integral to be finite it means that the integrand 
is peaked around momenta of order gioop ~ k. Using dimensional analysis we therefore conclude 
that 

PL-loops-finite ~ ^ ^^J l^J ' ^''^ 

in agreement with (42). Notice that this is a small correction only for k <^ /cnl. As discussed and 
used in [2], the EFTofLSS provides a manifestly convergent perturbative expansion in /c/Zcnl after 
renormalization has been performed. At this point each loop contributes a finite hierarchically- 
smaller correction in k/k^i^ -C 1. We also see that we are really interested in the finite contribu- 
tions. This means that the numerical precision in evaluating the integrals should be rather high. 
We now proceed to the numerical evaluation of the 2-loop integrals. 

3.2 Numerical Evaluation and Renormalization 

3.2.1 n = -|tilt 

Because of the UV divergence in the 2-loop integrals for the case n = —3/2, we will need to 
evaluate the 2-loop integrals with a cutoff A and cancel the A dependence with the speed of sound 
counterterm. This is given by 

p _ ^counter [ ^ ) / ^ ) p (/\f;\ 

-'2-loop counterterm /r> \t; 1 j I 1 j I -^ 11 ) V^^/ 

(27r)5 V«;nl/ yk^hj 

where the coefficient c^ounter is carefully chosen to cancel the UV divergence in P2-ioop- 

We will evaluate the integrals numerically with a finite A and with k <^ A;nl- We expect 



p 1 

-'2- Loop 



(27r) 



1/A I k\ f k \ „ ,, ,. „. . k 



+Ci -r -; — Pii + subleading finite terms in 

' A / V fcNL / A 
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Figure 1: Left: The 2-loop power spectrum P2-I00P for various values of the UV cutoff A and IR cutoff 
kmm- We show A = 2A;nl in blue, 10 /cnl in red, 20 /cnl in yellow, and SO/cnl in green, The result is /cmm 



independent, as we demonstrate in detail in Fig. 4, so here we plot only results with k„ 



5 X IQ- 



as 



the results with k^am = 5 x 10~^ are visually indistinguishable. Right: same as above but on a logarithmic 
scale. When the UV cutoff is high enough so that finite terms are negligible, we can see that all the curves 
have the same slope, corresponding to k^'"^ dependence, and that they depend linearly on A. This agrees 
with what is expected from the divergent term in (47). 



The factors of 27r count the phase-space suppression of higher loop contributions. In the last line, 
we have added terms that depend on A in a way that converges to zero as A — )■ 00. These terms 
are negligible only in the limit A — > 00. If one does the calculation at finite A, or one measures the 
coefficients of the EFTofLSS from A^-body simulations, these terms might need to be included, 
depending on the value of k/A. These terms are potentially important for us, as we evaluate the 
loops with a finite cutoff. 

If we evaluate the integrals with a very high A, it is then numerically difficult to have enough 
precision to extract the subleading finite term. We proceed therefore in the following way. We 
first evaluate the 2-loop integral with a very high-A, in such a way that we can neglect all the 
subleading terms and we can extract the coefficient c^ by measuring low-fc scaling behavior. In 
Fig. 1 we plot the value of the 2-loop power spectrum for various values of A and an IR cutoff fcmin- 
We see that for values of A high enough so that the finite terms are negligible, the result goes 
as Afc^/^, as it should, with no additional log(A) dependence, as anticipated by considerations 
related to the absence of counterterms in the EFTofLSS. This will be more clearly visible after 
we remove the linear divergence. We also note here that the result is k^i^ independent, as we 
demonstrate in detail in Fig. 4. 

After we have measured c^, we set 



^counter 



•^1 T^ '-"^counter 



X 



(4J 
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where ^Ccounter represents the finite (A- independent) contribution of the term in (46). As described 
above, we find cf by evaluating P2-I00P with a high cutoff A = 80 /cnl, and we fit the result with a 
k and A dependence of the form 



cf 



(27r) 



A 



NL 



k 



NL 



Pn 



(49) 



in the interval 0.04 /cnl < A; < 0.16 /cnl, so that the finite corrections, that become relatively 
smaller for A; — )■ 0, are under control. We take fcmin = 0.005 /cnl- As demonstrated later in Fig. 4, 
decreasing kmin by order of magnitude is completely negligible. We numerically find 



c^~ 



-2.277 



(50) 



which is order one as expected. We then evaluate the same P2-I00P with four lower cutoffs A = 
2 /cnl, 10 ^NL, 20 /cnl and 80 /cnl- After subtracting the divergent term, we fit the remaining part 
with the following functional form 



(27r)^ 



finite 

1 



+ C 



1/A 



+ cl^' 



A 



k 



NL 



2(3+n) 



Pu 



in the k interval 0.04 /cnl "^ k < A;nl- We numerically find 

-0.874 , c^^ ~ 9.005 , c^^^ ~ -5.745 , 



finite ^ 
1^1 



n 



-1.496 . 



(51) 



(52) 



The coefficients are of order one, as expected. Importantly, the measured scaling n agrees with 
the expected one n = —3/2. The fit comparison with the numerical data after the counterterm 
subtraction is shown in Fig. 2. Of course by including terms higher order in k/A one can improve 
the fit higher in k/kNL, but the above is sufficient to demonstrate the convergence of the IR- 
behavior of the integrand. 



3.2.2 



n 



I tilt 



For a tilt of n 



the UV situation is even more straightforward as there is no dominant 



dependence on the UV cutoff, the integral being UV convergent. This universe is a particularly 
stringent test of the IR behavior however. While it should be IR finite, there can and is measurable 
IR cutoff (fcmin) dependence, in a convergent way: 



1 1 



2-Loop, n=-5/2 



{2nfkl 



NL 




(53) 



This is entirely consistent with the idea of IR divergence occurring with n = —3. Consider for 
example the case of n = — 3 + e, with < e <^ 1. The dependence on kyam should be of order 



Finite P2-I00P for scaling universe Piioc(k/k^i^) ^^^ 




0.4 0.6 

k/^NL 



Figure 2: The points represents the finite contributions to P2-I00P for A = 2 A;nl (blue) A = 10 /cnl ('^^d), 
A = 20 /cnl (yellow) and A = 80 A;nl (green) obtained after the divergent part has been subtracted. As 
we go to low enough k 's, the three curves approach the same one, as indicated by the continuous fitting 
functions of Eq. (51), plotted as solid lines. As we move to higher k 's the A dependence, strongest for 
the A = 2 ^NL case, becomes clearly visible. 



{krain/kY , whlch is not a very subdominant correction for a given ratio femin/^ if e is small enough. 
To this point, under various IR cutoffs, we find the n = —5/2 tilt universe scales as follows: 



Dfit 
2-Loop, n=-5/2 



(271)5^3^ \knL 



k 



-3/2 



0.051-0.159 



/Cri 



1/2 



+ 0.114 



fCn 



(54) 



We present the data, and comparison to the fit in Fig. 3. The agreement is very good. 

4 Summary 

In the upcoming era where large scale structure surveys will likely be the leading source of cosmo- 
logical information, it is crucial to make reliable predictions for as many of the modes as possible 
in order to maximize the information that can be extracted. Since the phase space density of 
modes goes as k^, with k being the wavenumber of a mode, most of the information is contained 
in high wavenumbers, where non-linear corrections become important. In our current universe, 
there is a range of /c-modes, roughly between 0.1 hMpc~^ and 1 hMpc~^, where the non-linearities 
of dark matter clustering are large enough to make it necessary for them to be taken into account. 
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IR— cutoff dependence of Scaling Universe Pi i oc(k/iNL) 



n=-5/2 deviation: (data - fit)/ data 




Figure 3: Left: 2-loop power spectrum P2-I00P for ideal scaling of tilt n = —5/2 for various values of the 
IR cutoff krain- This univcrsc is UV cutoff (A) independent. Right: We demonstrate the functional form 
of the IR-dependence by comparing to the fit given in eq. (54)- 



but also small enough to be amenable to an analytical treatment. Given the phase space density 
of modes, this interval in wavenumber is actually the most important for observations. 

With the development of the Effective Field Theory of Large Scale Structures (EFTofLSS) [1, 
2], the theoretical framework where a manifestly convergent perturbative expansion in k/k-^i^ has 
been established, where k is the wavenumber of the mode of interest, and /cnl is the one associated 
with the non-linear scale. By computing higher order corrections, adding suitable counterterms 
order by order, and renormalizing the theory, more accurate predictions in k/k^i^ can be made, 
approaching the non-linear scale /cnl- 

Given how much phase space is focussed at the UV and the absence of formerly reliable 
analytical techniques, the actual value of /cnl in our universe at order-one is still unclear. Having 
/cnl = 0.5 /iMpc"^ or 1 hMpc"^ makes a very large difference in the amount of available data we 
can use for cosmology. 

We suggest /cnl should be the scale at which the EFTofLSS stop converging. To identify this, 
higher order calculations beyond the ones first studied in [2] need to be computed. The next step 
is to compute the 2-loop correction. Beyond the standard numerical challenges, one is faced with 
two additional subtleties. First, each of the standard two-loop diagrams has several IR divergences 
in different soft regions of the integrand that cancel only when all the diagrams are summed over. 
This makes the numerical evaluation much harder, as the largest contribution to each diagram 
is actually unphysical and will need to cancel out from the final answer. Second, the leading 
contribution of the IR-independent part is degenerate with the speed-of-sound counterterm in the 
EFTofLSS, and so is uncalculable. For this reason, the physically relevant 2-loop term is the 
next-to-leading IR-finite term. 

In this paper we have presented a procedure which re-arranges all the four standard 2-loop 
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diagrams into only one integral, and we manipulate the integrand in such a way that the only 
possible IR divergences appear when either both the internal momenta or only one of the two goes 
to zero. Since IR divergences are guaranteed to cancel from the final answer, our procedure enforces 
the cancellation at the level of the integrand. Now this 'IR-safe global integrand' is IR integrable 
in any soft limit, and this happens both for the leading IR divergences where both the internal 
momenta go to zero, and also for the sub-leading ones where only one of the internal momenta go 
to zero. This makes the numerical integration much easier, as the integral is automatically of the 
size of the IR- independent part, without there being cancellation between regions of the integrand 
far-away in phase space. 

After implementing our IR-safe integrand we have computed the 2-loop power spectrum in 
scaling universes with a simple linear power spectrum of the form k"'. They have the property of 
being particularly simple, so that the form of the final answer can be derived using dimensional 
analysis. At the same time these form the most strident tests from the IR-point of view, as the 
standard integrand does not go nicely to zero in the soft limit. In particular, we have explored first 
the case n = —3/2, which has the property of being UV divergent. We evaluate the integrals with 
a UV cutoff A, we renormalize the leading divergence, and we investigate the residual dependence 
on the IR cutoff /Cmin and verify final answer agrees with the behavior we predicted from dimen- 
sional grounds. Satisfactorily, we find that the dependence on fcmin is subpercent, and that the 
residual analytic dependence of the 2-loop power spectra after renormalization perfectly matches 
our expectations. Subsequently we have tested the case n = —5/2, which does not have UV 
divergences, so that the result of the integration is automatically physical and no renormalization 
is needed, but that has a potentially larger IR divergence. In this case as well, we have found that 
the answer depends correctly on k^in, with corrections proportional to (k^m/kY^'^, which goes 
to zero as fcmin — > 0. We have also tested our calculation against different numerical integration 
techniques (Montecarlo and Quasi-Montecarlo), finding that the result does not depend relevantly 
on the numerical technique used. 

We have compared the performance of our technique for computing the 2-loop power spectrum 
with the standard techniques used in some publicly available codes. We have found that without 
the implementation of our IR-safe integrand, these code are unable to reliably compute the 2-loop 
power spectra for scaling universe, and probably also for our current universe. 

In a forthcoming paper [11], we will apply this integration technique to give the 2.5-loop 
prediction of the power spectrum from the EFTofLSS, where we will see that the EFTofLSS 
predicts the power spectrum with percent precision up to roughly k ~ 0.5 hMpc~^. 
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Appendix 

A Numerical Checks 

A.l Independence of IR cutoff and of method of integration 

In this appendix we provide some details of the implementation and some explicit checks of 
numerical convergence. 

All numerical integrations of the global IR-safe integrand were performed using the CUBA li- 
brary [13] (v. 3.0). The implementation of our IR-safe integrand used the recursive implementation 
of the kernels given in the Copter library [4], but combined as given above in eq. (33). 

First we demonstrate that our results are independent of the IR cutoff. In Fig. 4, we plot 
the relative difference for P2-I00P obtained for different values of krain- We see that over the range 
explored, the /cmin-dependence is at the less than percent level, becoming of order 0.3% at fc ~ /cnl- 

We demonstrate our results are independent of the method of integration in Fig. 5. We find 
that there is sub-percent relative difference in P2-I00P computed for A = 2, using two methods, 
"quasi-Monte- Carlo" and "random Monte-Carlo," for the Vegas algorithm in CUBA. Both were 
set to terminate when achieving either a relative (estimated) accuracy goal of 10"'^, or reaching a 
maximum of 10^ integrand evaluations. 

A. 2 Comparison with available 2-loop results 

In this section we compare our results with other publicly available numerical tools and their 
associated 2-loop integrands. We do this to put into evidence the importance of performing the 
numerical integration with the IR-safe integral ^ . This is relevant not only because it is difficult to 
cancel the IR divergences numerically, but also because, as for the case of n = —3/2 and for the 
ACDM universe, we are interested in the subleading part of the integrand, the leading part being 
degenerate with a counterterm of the EFTofLSS that is not calculable within the EFT anyway. 



^We stress the codes we will use here for comparison are not designed for scaling universes. It might well be 
that some of these IR-issues are numerically irrelevant for these codes in the case of the true universe. 
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Figure 4: Fractional difference between P2-ioop,n=-3/2 f(^f ^min = 5 x 10 ^ /cnl and /cmin 
for A = 10. We see that for k <1, the dependence is less than 1/3 of a percent. 
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We first consider the Copter evaluation available in Ref. [4]. This approach simply uses a 
quasi- Monte Carlo method method to separately integrate the four diagrams P51, P42, P33 and 
P33 . The comparison is unfortunately straightforward, as we were unable to get Copter to return 
anything but "NAN" error, which means "not-a-number" error, for scaling universe input, even 
with finite UV and IR cutoffs. While there may be a simple fix to get numbers out for comparison, 
we did not find one. 

We then compare with the RegPT approach of Ref. [10]. While this is built to implement a 
regularised perturbation theory, one of the available outputs is the 2-loop SPT power spectrum, 
which is the one we are interested in for the purpose of comparison. This code does not treat the 
IR divergences in the way we do, but it is still able to complete calculations for scaling universes.^ 

In Fig. 6, we plot the P2-I00P from RegPT with fcmm = 0.005 A;nl (solid points) and fcmin = 
0.0005 /cnl (empty points) for the same values of A (2fcNL, 10 /cnl, 20 /cnl, and 80 /cnl) as in Fig. 2. 
With the larger IR cutoff, we observe the correct Ak^'"^ behavior, but begin to see scatter of the 
calculated points around the smooth fit lines at higher k. For the smaller IR cutoff, the compu- 
tation becomes unstable. Even for the larger IR cutoff, numerical noise becomes overwhelming 
when we cancel the UV-divergent piece with the appropriate counterterm and plot the leftover 



^One should note that as of the time of writing, RegPT requires an earher version of CUBA: (v. 1.5). We verified 
that the IR-safe global integrand results are unaffected by differences in CUBA versions. 
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Figure 5: Plot of the deviation between the "quasi- Monte- Carlo" and "random Monte- Carlo" results 
from CUBA'S Vegas integration algorithm, for the IR safe integrand for scaling universe n = —3/2. We 
see that it is sub-percent for the range considered. 



finite piece. The code's default precision is apparently insufficient to provide meaningful results 
for UV- convergent part of the P2-ioop; as evidenced by the mismatch between the numerical results 
from RegPT and the fitting functions obtained from our IR-safe integrand (plotted as solid lines 
in the right-hand frame). 

We reiterate that the physical information of P2-I00P is contained in the UV-convergent part. 
Therefore any numerical calculation of P2-I00P for scaling universes should not only reproduce the 
A-dependence seen in the left-hand frame of Fig. 6 (or the corresponding version for different 
tilts), but should also be able to resolve the behavior of the UV-convergent piece with the wished 
precision. 

It is also possible to increase the precision of the integration methods used by RegPT and 
examine the results. We show such a computation in Fig. 7, where the code's runtime has increased 
by a factor of ~100 as a result of the heightened precision. While the results with larger k^[^ 
are qualitatively consistent with those of our IR-safe integrand (although there are noticeable 
deviations for A = 10 ^nl and 20/cnl), there is still strong IR-dependence that creates numerical 
issues even at higher than default precision. These issues become even more severe as fcmin is taken 
lower and lower. In principle, running RegPT using arbitrarily high numerical precision would 
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Figure 6: Left: P2-I00P output of RegPT for a scaling universe with n = —3/2, using k^^in = 0.005A;nl 
(solid points) and /cmin = 0.0005 /cnl (empty points). We show A = 2 /cnl (blue), IO/cnl (fed), 20A;nl 
(yellow), and 80 /cnl (green) and omit two small-kinm curves to avoid clutter. The solid lines are fits to 
the dominant A k^'"^ behavior of each set of points — notice the scatter of the points about the smooth fits 
at higher k, as well as the strong k^\n- dependence, signalling that IR-divergences are not properly treated. 
Right: P2-I00P after having subtracted the UV-divergent piece of each curve. The dotted lines connect the 
computed points as a visual aid, while the solid lines are the fitting functions from Eq. (51). Numerical 
instabilities prevent RegPT from reliably indicating the form of the terms (compare with Fig. 2). 



likely allow the user to recover the detailed forms of the finite part of P2-ioop- However, a more 
robust approach would implement the calculations in a way that makes IR divergences manifestly 
absent before the integrals are evaluated, as with our IR-safe integrand. 
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